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ABSTRACT 

According to ray theory, regions exist in an upward refracting atmosphere where no sound should be 
present. Experiments show, however, that appreciable sound levels penetrate these so-called shadow 
zones. Two mechanisms contribute to sound in the shadow zone: diffraction and turbulent scattering of 
sound. Diffractive effects can be pronounced at lower frequencies but are small at high frequencies. In 
the short wavelength limit, then, scattering due to turbulence should be the predominant mechanism 
involved in producing the sound levels measured in shadow zones. No existing analytical method 
includes turbulence effects in the prediction of sound pressure levels in upward refractive shadow zones. 

In order to obtain quantitative average sound pressure level predictions, a numerical simulation of the 
effect of atmospheric turbulence on sound propagation is performed. The simulation is based on scattering 
from randomly distributed scattering centers ("turbules"). Sound pressure levels are computed for many 
realizations of a turbulent atmosphere. Predictions from the numerical simulation are compared with 
existing theories and experimental data. 




INTRODUCTION 

Solar heating of the ground produces strong temperature gradients in the air just above the surface of 
the Earth. Since the speed of sound is proportional to the square root of the temperature, sound will 
follow upwardly curved paths in every direction from a source. The stronger the temperature gradients 
involved, the shorter the distance to what is properly called a shadow zone, since no direct or reflected 
rays can penetrate into this region. Figure 1 is an illustration of an upward refractive shadow zone where 
the source is at a height h s and the radius of curvature of the limiting ray is The edge of the shadow 
zone is delineated by the so-called limiting ray which grazes the ground. 

In a similar fashion, sound traveling upwind is curved upwards due to the strong wind gradients 
near the ground and a shadow zone is also formed. In the case of wind, however, the effect is not 
isotropic because of the vector nature of the wind velocity and the rays are actually bent downward for the 
sound propagating downwind. 

Two mechanisms contribute to the magnitude of the sound levels measured in shadow zones: 
diffraction and the turbulent scattering of sound. Pierce 1 describes the solution for a linear sound speed 


245 


gradient in terms of residue series for the pressure in the shadow zone. He examines the cases of a hard 
boundary and a pressure-release surface and gives approximate solutions when the source and receiver are 
above the creeping wave layer height, defined as (Rg/2 where R^. is the radius of curvature of the 

limiting ray and k 0 is the wavenumber. 

Daigle et al. 2 made use of the above two approximate solutions in an effort to fit the data they 
collected over an asphalt airport runway and over a grass-covered strip near the runway, the latter 
approximating a pressure-release surface at frequencies greater than 500 Hz and the former approximating 
a hard boundary. They found that the hard boundary data was well explained by Pierce's approximate 
solution for that case and that the data up to 1000 Hz over the grass-covered ground was satisfactorily 
explained by Pierce's approximate solution for a pressure-release surface. 

The approximate solution leads to large errors in the effective source levels for sources close to a 
pressure-release or finite impedance ground, as was the case with Daigle's data. A complete discussion of 
this problem can be found in the paper by Raspet and Franke. 3 In a later paper. Berry and Daigle 4 used 
the complete residue series solution and again compared the above data. They found that the data at 250 
Hz still agreed well with the predictions of diffraction theory. But the data was well under-predicted by 
diffraction theory at 500 Hz and especially at 1000 Hz. The predictions from the full residue series 
solution are shown in Fig. 2, which is a reproduction of Fig. 13b of Ref. 4. 

The role played by atmospheric turbulence in the insonification of shadow zones has escaped 
analytical formulation. In an effort to obtain a quantitative estimate of the extent to which atmospheric 
turbulence raises the sound levels in a shadow zone, Gilbert et al. 3 used a parabolic equation method to 
numerically simulate sound propagation in a turbulent atmosphere. They compared their predictions for 
upward refracting conditions with experimental results of Wiener and Keast.6 The numerical predictions 
involved the calculation of the sound pressure magnitude for a particular realization or "snapshot" of 
turbulence, while the results of Wiener and Keast were expressed in terms of average sound pressure 
levels. Nevertheless, Gilbert et al. were able to duplicate the apparent range independence of excess 
attenuation characteristic of the experimental data at ranges as great as 1 km. 

In this paper, we present the average sound pressure levels in an upward refractive shadow zone 
predicted by a scattering center based numerical simulation. The main features of the numerical solution 
are reviewed and the modifications necessary to adapt it to an upward refractive atmosphere are discussed. 
Sound levels are computed for over 500 realizations of the turbulent atmosphere. Predictions from the 
numerical simulation are then compared with experimental data taken by Daigle et al. 2 

MAIN FEATURES OF THE NUMERICAL SIMULATION 

Although the details of the numerical simulation were given in an earlier paper* the main features are 
repeated here so that the reader may have a better idea of the type of calculations involved. Following the 
model of de Wolf, 7 we construct an ensemble of isotropic, irrotational scattering centers which we call 
"turbules." If |i is defined as the change from unity of the index of refraction, a given turbule is assigned 
the refractive profile 

* Walton E. McBride, Henry E. Bass, Richard Raspet, Kenneth E. Gilbert, "Scattering of sound by 
atmospheric turbulence," submitted to J. Acoust. Soc. Am., Feb. 1990. 
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( 1 ) 


[L(r,s) = qie' r/S 

where < 7 j is the value of |1 at the center of the spherically symmetric turbule and s is the 1/e contour of 
the scattering center and can be considered to be its effective size. The value of < 7 , and the probability 
distribution of turbule sizes depend in general on the particular functional form chosen for the correlation 
function of the fluctuations of the index of refraction. If the correlation function is chosen to have the 
Gaussian form, 


) = ^ e ' r/L > (2) 

where <p 2 > is the variance and L is the correlation length, then the size spectrum is a delta function 
implying that all the turbules have the same size. 



The value of < 7 ,- for this particular form of the correlation function is given by 


(3) 



(4) 


and is inversely proportional to the turbule number density p^. An upper value of p^j of about half the 
overlap density is necessary so that the turbules will be separate entities. With single scattering, sound 
scattered from a particular turbule reaches the receiver downfield with negligible scattering by other 
turbules located between that particular turbule and the receiver. From Eq. (4) the product < 7 ,- 2 p N is a 
constant whose value depends on the independently measured micrometeorological variables <p 2 > and L. 
There is, therefore, a certain latitude in the value of p^. Decreasing p N will result in a greater value for 
l< 7 ,l. Although a lesser number of turbules result from a decrease of p N , the predictions of the numerical 
simulation are statistically steady as the turbule number density is decreased from an upper limit of half the 
overlap density. 

Initially in the development of the simulation, the first Bom approximation to scattering was used to 
determine the scattering effect of each turbule. In practice, the evaluation of the scattering integral is 
performed by assuming that both source and receiver are far away from the scattering region; thus first 
order terms in the phase are sufficient. Because some of the turbules in the numerical simulation are close 
to the source or receiver, second order terms in the phase were kept in the scattering integral. The total 
pressure at a receiver downfield is, then, the sum of the direct and scattered spherical waves. For only one 
turbule in free space, this is: 
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(5) 


p(R) = 


R 



j^2 s 3 e «'*(*> n-o,) 


1 ^ e -CJfcV/4 


where 


C = (l-cos ©o) 2 + sin 2 ©o 
and 



( 6 ) 


(7) 


In the above equation, k is the wavenumber, s is the effective size of the turbule, © 0 is the angle 
between the incident and scattered directions, R is the distance between source and receiver, while r st is the 
distance between source and turbule center, and r^ is the distance between turbule center and receiver. 

Note that the usual Bom scattering term is recovered when a = 0 and the first term in Eq. (6) is 
dropped. Even with this improved evaluation of the Bom scattering integral, the distance from turbule to 
source or receiver cannot be less than about twice the radius of the turbule. Consequently, "buffers" of a 
turbule's diameter were placed in front of the source and receiver where no turbules were allowed. 

The numerical simulation using the first Bom approximation to scattering was then compared to 
theoretical expressions due to Karavainikov^ for the log-amplitude and phase variances of the pressure 
fluctuations. It was found that good agreement was reached whenever the wave parameter D (=R/&L 2 ) 
was greater than 1. As shown in Fig. 3, the log-amplitude variances as predicted by Karavainikov are 
independent of frequency when D < 1, a result also obtained by Bergmann using geometrical optics. 

In an effort to reach better agreement in the geometrical optics region, the Rytov approximation used 
by Karavainikov was incorporated into the numerical simulation. The Rytov method consists of 
approximating the field at the receiver by 


p R (f) = p 0 (r)eV R ®; 


( 8 ) 


whereas the Bom approximation is written: 


p B (r) = p 0 (f) + ? B (?X 

There is a simple relationship between the first Rytov and first Bom approximations: 


(9) 


Po(?) 


(10) 


and it was, therefore, a simple matter to incorporate the first Rytov approximation into the numerical 
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simulation. The results are shown in Fig. 4 and good agreement is obtained throughout the range of the 
wave parameter D. 

As can be seen from the above comparisons with Karavainikov's analytic curves, the first Rytov 
approximation is superior to the first Bom approximation for an unbounded medium without refraction. 
When refractive conditions are introduced, however, the formation of shadow zones becomes possible. In 
the shadow zones, p 0 is 0 and the first Rytov approximation cannot be used. Recourse must be made to 
the first Bom approximation and the wave parameter D must be greater than 1 for the numerical simulation 
to be valid in accordance with the results of Fig. 3. 

The next step is the inclusion of the ground. An immediate consequence of the existence of a 
boundary is the presence of three additional paths by which sound can propagate to the receiver. There 
now exist four single scatter paths that connect the source and receiver: 

1 . source-turbule-receiver, 

2. source-turbule-ground-receiver, 

3 . source-ground-turbule-receiver, 

4 . source-ground-turbule-ground-receiver. 

The last three paths all interact with the ground and, therefore, a model of the effect of the ground on the 
sound wave was also included in the numerical simulation. 

The algorithm proceeds as follows. Values of <(i 2 > and L are given from independent micro- 
meteorological measurements. From these, the value of s and are obtained using the above equations. 

A scattering space, which will enclose thousands of turbules, is defined with buffers in front of the source 
and receivers of widths equal to about the diameter of a turbule. The turbules are assigned positive or 
negative signs for their value of q v The sound pressure at the receiver is calculated for this particular 
arrangement of turbules, and the result is referred to as a realization. Then each turbule is given random, 
small increments in its Cartesian coordinates. The sound pressure at the receiver is recalculated, resulting 
in another realization. The process is repeated for as many realizations as are necessary for the statistics to 
stabilize. We have found that 500 realizations are sufficient. Average sound pressure levels can then be 
obtained from the 500 stored values of the sound pressure. It should be mentioned that any other desired 
statistical quantity can be obtained, such as structure and correlation functions, as well as the variances of 
the log-amplitude and phase fluctuations. 

The inclusion of a sound speed gradient requires two modifications: the rays are now curved and the 
value of the wavenumber k is no longer constant along a ray. In order to obtain a closed form solution 
for the equation describing the rays, a linear sound speed gradient was assumed. As is well known, a 
consequence of this assumption is that the ray paths are arcs of circles. 

Because of the curvature of the ray paths, each path must be tested to see whether the source and 
each turbule can be joined together, as well as each turbule and the receiver. If either segment of the total 
path cannot be linked, that particular turbule's contribution is discarded. It was found that about 15% of 
the turbules were eliminated in this way for the experimental data to be described later. 

The last correction necessary is the calculation of an effective wavenumber k e for each ray path. 

This required the computation of the length of each path, as well as the travel time along that path. Their 
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ratio gave an effective sound speed c e along that path, and the effective wavenumber was then given by 
oVc e . 

COMPARISON TO DATA 

In order to implement the numerical simulation, the statistical properties <|i 2 > and L of the 
turbulent atmosphere and the impedance of the ground are required. The former was given in the article by 
Daigle et al. 2 as <p 2 > = 6 x 10' 6 and L = 1.6 m. The impedance had to be approximated because the 
article mentioned above did not specify a particular impedance model. To estimate the impedance, a 
residue series solution developed by Raspet and Franke^ was used to match the curves of Fig. 2 as closely 
as possible at all three frequencies. The particular impedance model used was a four parameter model 
developed by Attenborough.^ A shape factor ri of .750, a shape factor ratio Sj - of .875, a porosity Q of 

.675, and a flow resistivity a of 330 cgs rayls give the results shown in Fig. 5. 

The numerical simulation was performed with the above parameter values for 500 realizations. Rms 
sound pressure values were computed and divided by the pressure doubling plus attenuation factor 
(2 e'QK/R) as was done in Daigle's presentation of his experimental data. The source was given a height 
of 0 m, and six receiver positions were used 230 m downrange at heights of 0.25 m, 0.50 m, 1.0 m, 2.0 
m, 4.0 m, and 7.0 m in order to sample the vertical behavior of the sound pressure levels in the shadow 
zone. Because we are forced to use the first Bom approximation when the receiver is in a shadow zone 
(pQ = 0), the possible values of the wave parameter D must be checked to see that they will be greater than 
1. Since D = R/fcL 2 with R = 230 m and L = 1 .6 m, it is sufficient to check the value of D for the greatest 
frequency of interest. With 1000 Hz, k ~ 18.4 and, therefore, D = 5 which is well above the minimum 
value of 1. The sound pressure levels predicted by the numerical simulation (expressed in dB) are 
compared to Daigle's experimental data in Fig. 6 for the three frequencies involved. Notice the similarity 
in shape which the curves representing the turbulent scattering contribution share with the curves that are 
typical of diffraction theory predictions. As can be seen by comparing Figs. 2 and 6, it appears that for 
this experiment the contributions from turbulent scattering and diffractive effects are about equal at 250 
Hz. However, turbulent scattering becomes the major contributor at 500 Hz. At 1000 Hz, turbulent 
scattering is the predominant mechanism behind the increased sound pressure levels measured in the 
shadow zone. 


CONCLUSION 

Quantitative predictions for the average sound pressure levels in a refractive shadow zone have been 
presented. Good agreement was reached with experimental data collected by Daigle et al 2 in a shadow 
zone caused by temperature gradients. It seems that the use of a simple linear sound speed gradient is a 
good approximation to the real sound speed profile directly above the ground for the moderate ranges 
involved in this study. At the longer ranges investigated by Gilbert et al. it was necessary to use a 
logarithmic sound velocity profile to obtain accurate predictions.^ The relative contributions of diffraction 
and turbulent scattering have been examined and graphically displayed. The dominant mechanism which 
dictates sound levels in shadow zones at higher frequencies is scattering due to turbulence. 
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Fig. 1 Shadow zone formation for a sound speed profile that decreases linearly with height. 



Fig. 2 Comparison of measured sound levels (points) with predictions based upon diffraction into the 
shadow zone, taken from Ref. 4. Solid circles are for 1000 Hz, triangles for 500 Hz, and 
diamonds for 250 Hz. 


252 





10 


o 

§ 

1.0 - 

*s 

> 

o 

cn 

10 1 - 

a 

-C 

O* 

10 2 - 

a 

'E, 

10-3- 

e 

< 

bO 

UH- 

3 

10’ 5 - 



0.1 0.4 0.7 1.0 4 7 10 


Frequency (khz) 


Fig. 3 Comparison of the numerical simulation using the first Bom approximation (O for phase variances; 

X for log-amplitude variances) with Karavainikov's analytical results ( for phase variances; 

for log-amplitude variances). 



Fig. 4 Comparison of the numerical simulation using the first Rytov approximation (o for phase 

variances; X for log-amplitude variances) with Karavainikov's analytical results ( for phase 

variances; for log-amplitude variances). 
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Fig. 5 Comparison of data with the prediction of diffraction theory using the residue series solution. Data 
points are described for Fig. 2. 



Relative dB Level 


Fig. 6 Comparison of the prediction from the numerical simulation with experimental data taken in 
shadow zone. Experimental data is the same as in Fig. 5. 
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